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From  B-Spline  Mesh  Generation  to  Effective 
Visualizations  for  the  NIST  Digital  Library  of 
Mathematical  Functions 

Bonita  Saunders,  Qiming  Wang 

Abstract.  We  discuss  the  use  of  tensor  product  B-splines  to  generate 
grids,  or  meshes,  to  facilitate  the  plotting  of  high  level  mathemati- 
cal functions  for  the  National  Institute  of  Standards  and  Technology 
(NIST)  Digital  Library  of  Mathematical  Functions  Project.  The  plot 
data  is  placed  inside  a web  based  format  such  as  VRML  (Virtual 
Reality  Modeling  Language)  to  create  interactive  visualizations  that 
allow  users  to  carefully  examine  complicated  function  features  such  as 
zeros,  branch  cuts,  poles  and  other  singularities.  We  discuss  the  grid 
generation  technique  and  look  at  its  effectiveness  in  creating  accurate 
visualizations. 


§1.  Introduction 

The  National  Institute  of  Standards  and  Technology  (NIST),  a U.S.  gov- 
ernment agency  under  the  Department  of  Commerce,  is  in  the  midst  of  de- 
veloping the  NIST  Digital  Library  of  Mathematical  Functions  (DLMF),  a 
web-based  collection  of  high  level,  or  special,  mathematical  functions  that 
will  replace  the  widely  used,  but  outdated,  National  Bureau  of  Standards 
Handbook  of  Mathematical  Functions  published  in  1964  [1],  The  website 
will  include  formulas,  methods  of  computation,  references,  and  links  to 
software  for  over  forty  functions.  It  will  feature  interactive  navigation,  a 
mathematical  equation  search,  2D  graphics,  and  dynamic  interactive  3D 
visualizations. 

This  paper  discusses  our  use  of  a mesh  generation  technique  to  facil- 
itate the  plotting  of  functions  for  the  3D  visualizations.  Often,  the  com- 
plicated nature  of  a special  function  is  clearly  indicated  by  its  domain, 
which  may  be  irregular,  discontinuous,  or  multiply  connected.  By  modify- 
ing an  algebraic  tensor  product  spline  mesh  generation  algorithm  that  we 
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originally  designed  for  problems  in  aerodynamics  and  solidification  theory, 
we  were  able  to  create  boundary/contour  fitted  computational  grids  that 
allowed  us  to  capture  key  function  features  such  as  zeros,  poles,  branch 
cuts  and  other  singularities. 

Although  commercial  packages  often  have  many  built-in  special  func- 
tions, their  3D  plots  are  usually  over  a rectangular  Cartesian  mesh,  lead- 
ing to  poor  and  misleading  graphs.  Also,  when  function  values  lie  outside 
the  range  of  interest,  many  packages  have  trouble  properly  clipping  the 
function  surface  to  produce  an  accurate  graphical  representation.  Further- 
more, even  when  the  plot  looks  satisfactory  inside  a package,  it  may  be 
completely  unacceptable  when  the  data  is  transformed  to  a format  more 
suitable  for  display  on  the  web. 

We  examine  these  problems  and  others  that  arise  when  trying  to  dis- 
play dynamic  3D  graphs  of  special  functions  on  the  web.  We  show  how 
our  mesh  generation  algorithm  eliminates  or  lessens  the  severity  of  many 
problems,  talk  about  progress  in  making  the  method  more  adaptive  and 
discuss  other  possibilities  for  improvements. 

§2.  Constructing  3D  Graphics  for  a Digital  Library 

The  NIST  DLMF  will  consist  of  approximately  forty  chapters  authored 
by  special  function  experts  throughout  the  U.  S.  and  abroad.  The  number 
and  location  of  graphics  for  each  chapter  is  determined  by  consultation 
with  the  authors  and  DLMF  editors.  The  accuracy  of  the  plotting  data 
is  being  verified  by  plotting  the  function  with  at  least  two  different  meth- 
ods, using  commercial  packages,  publicly  available  codes,  or  the  author’s 
personal  codes.  A key  concern,  however,  is  whether  or  not  the  displayed 
plot  accurately  represents  the  graph  of  the  function. 

Commercial  packages  typically  render  2D  plots  very  well.  Discontinu- 
ities are  handled  automatically  or  fairly  easily  with  special  options.  If  a 
user  wants  to  restrict  the  vertical  range  of  the  function,  the  package  prop- 
erly cuts,  or  clips,  the  curve  so  that  only  points  within  the  desired  range 
appear.  This  is  often  not  the  case  in  3D.  Figure  1 illustrates  some  of  the 
problems  we  have  encountered  when  trying  to  use  commercial  packages. 
The  first  figure  shows  a plot  of  l/|r(z)|,  where  z — x + iy,  rendered  using 
a popular  commercial  package.  The  user  has  requested  that  the  function 
only  be  plotted  for  values  less  than  or  equal  to  6.  The  package  produces 
a flat  area,  or  shelf,  where  values  greater  than  6 are  set  to  6.  Although 
the  shelf-like  area,  which  has  nothing  to  do  with  the  function,  may  not 
concern  many  users,  it  might  be  confusing  to  students  and  others  unfamil- 
iar with  the  function.  The  user  can  request  that  the  shelf  not  be  shown, 
but  this  produces  a saw-tooth  area  which  can  be  even  more  misleading. 
Although  seasoned  users  can  use  alternative  commands  to  successfully  clip 
the  function  properly,  we  found  that  this  was  still  not  sufficient  for  our 
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requirements.  The  problem  is  that  the  computation  of  the  function  is  still 
over  a rectangular  Cartesian  mesh.  The  figure  looks  fine  when  viewed 
inside  the  package,  but  when  we  put  the  data  into  our  web  format,  the 
rectangular  mesh  cells  produce  a very  irregular  color  map.  For  our  web 
visualizations,  we  use  the  Virtual  Reality  Modeling  Language  (VRML),  a 
standard  3D  file  format  for  creating  interactive  web-based  visualizations 
[6],  The  second  figure,  rendered  using  the  same  package,  shows  the  modu- 
lus of  the  complex  gamma  function,  |r(z)|,  z = x + iy.  Both  figures  can  be 
improved  by  using  a much  larger  number  of  mesh  points,  but  large  data 
files  hamper  the  interactive  features  of  the  web-based  VRML  visualiza- 
tions. The  rendering  problems  can  be  eliminated  or  decreased  in  severity 
by  computing  the  function  over  a specially  designed  boundary /contour  fit- 
ted mesh.  The  next  section  discusses  the  tensor  product  B-spline  mapping 
we  have  used  to  produce  such  meshes. 


Fig.  1.  Potential  problems  such  as  bad  clipping  and  poor  resolution  of  poles. 


§3.  Grid  Generation  Mapping 

For  our  application,  the  difficulty  of  the  mesh  generation  problem  depends 
on  the  complexity  of  the  computational  domain  for  the  special  function. 
The  domains  range  from  simple  rectangles  to  complicated  multiply  con- 
nected domains  with  branch  cuts.  Our  primary  mesh  generation  technique 
is  based  on  an  algorithm  developed  by  one  of  the  authors  for  meshes  to 
be  used  in  solving  partial  differential  equations  (pdes)  related  to  areody- 
namics  and  solidfification  theory  [4],  [5].  The  algorithm  uses  an  algebraic 
grid  generation  system  defined  by  a mapping  T from  the  unit  square  I2 
to  a physical  domain  of  arbitrary  shape.  Specifically,  we  let 

l x(£,ri)  \ (£>??)  \ 

V y(^v)  ) V E™iE"=i )' 
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where  0 < £,  r/  < 1 and  each  BtJ  is  the  tensor  product  of  cubic  B-splines. 
Hence,  Btj(^,rj)  = Bl(£)Bj(r])  where  Bl  and  Bj  are  elements  of  cubic 
B-spline  sequences  associated  with  finite  nondecreasing  knot  sequences, 
say,  {sj}jn+4  and  respectively.  For  simple  boundaries  the  coeffi- 

cients atj  and  Bij  can  be  chosen  to  produce  a very  good  approximation 
to  a common  algebraic  generation  technique,  transfinite  blending  function 
interpolation.  In  short,  this  means  the  coefficients  are  obtained  by  evalu- 
ating T at  average  knot  values  as  discussed  in  [2] . This  shape  preserving 
approximation  reproduces  straight  lines  and  preserves  convexity  [2]  . For 
more  complicated  or  highly  nonconvex  boundaries  we  can  turn  the  tech- 
nique into  a variational  method  with  the  coefficients  chosen  to  minimize 
the  functional 


where  T denotes  the  grid  generation  mapping,  J is  the  Jacobian  of  the 
mapping,  W\ , W2  and  W3  are  weight  constants,  and  u represents  external 
criteria  for  adapting  the  grid  . Like  the  variational  method  of  Brackbill 
and  Saltzman  [3],  the  integral  controls  mesh  smoothness,  orthogonality, 
and  depending  on  the  definition  of  u,  the  adaptive  concentration  of  the 
grid  lines.  Hence,  when  solving  pdes,  u might  be  the  gradient  of  the 
evolving  solution  or  an  approximation  of  truncation  error.  Ideally,  for  our 
problem,  we  want  u to  be  based  on  curvature  and  gradient  information 
related  to  the  function  surface.  To  avoid  solving  the  Euler  equations  for 
the  variational  problem,  this  functional  is  approximated  in  the  computer 
code  by  the  sum 
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where  Jtj  is  the  Jacobian  value,  Uy  is  the  value  of  u,  and  Dottj  is  the  dot 
product  of  DT/d£  and  dT/drj  at  mesh  point  (£*,77 j)  on  the  unit  square. 
When  W3  = 0,  G is  actually  a fourth  degree  polynomial  in  each  spline 
coefficient  so  the  minimum  can  be  found  by  using  a cyclic  coordinate 
descent  technique  which  sequentially  finds  the  minimum  with  respect  to 
each  coefficient.  This  technique  allows  the  minimization  routine  to  take 
advantage  of  the  small  support  of  B-splines  when  evaluating  the  sums  that 
comprise  G. 
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An  application  of  the  grid  generation  algorithm  is  shown  in  Figure  2 
for  a puzzle  shaped  domain.  The  initial  grid,  constructed  using  linear 
Lagrange  polynomials  for  the  blending  functions,  is  shown  on  the  left. 
Note  that  the  grid  lines  overlap  the  nonconvex  boundary.  The  grid  on  the 
right  shows  the  mesh  obtained  after  the  spline  coefficients  are  modified  to 
minimize  G.  The  overlapping  grid  lines  have  been  pulled  into  the  interior. 


Fig.  2.  Initial  and  optimized  puzzle  grids. 


§4.  Results 

We  have  used  boundary /contour  fitted  grid  generation  to  produce  compu- 
tational grids  for  over  one  hundred  3D  visualizations  for  the  NIST  DLMF. 
The  plots  in  Figure  3 show  the  contour  boundary  and  mesh  for  the  func- 
tion l/|r(z)|.  To  obtain  the  boundary  we  found  the  contour  curves  for 
l/|r(z)|  = 6.  We  then  connected  the  curves  to  parts  of  a rectangle  to 
form  a closed  boundary  and  generated  the  boundary/contour  fitted  grid 
shown  on  the  right.  The  picture  on  the  left  in  Figure  4 shows  the  sur- 
face obtained  after  computing  l/|F(z)|  over  a purely  rectangular  grid  with 
values  outside  the  desired  range  set  to  6.  The  picture  on  the  right  shows 
the  result  when  we  compute  over  our  boundary/contour  fitted  grid.  In 
both  pictures  we  show  a snapshot  of  the  surface  rendering  obtained  after 
the  data  is  translated  into  VRML  format  for  viewing  on  the  web.  The 
boundary /contour  fitted  grid  properly  clips  the  function  and  improves  the 
smoothness  of  the  color  map.  Figure  5 shows  a grid  and  surface  for 
the  modulus  of  the  complex  gamma  function  |r(z)|.  The  top  and  bottom 
halves  were  generated  separately,  with  an  exponential  function  used  to 
concentrate  the  grid  points  near  y = 0.  The  surface  shown  on  the  right 
illustrates  how  precisely  the  contour  curves  satisfying  |r(z)|  = 6 were  com- 
puted. Once  the  coefficients  are  defined,  we  have  the  option  of  creating 
a coarser,  or  finer,  grid  simply  be  evaluating  fewer,  or  more,  points  on 
the  unit  square.  Also,  notice  that  the  grid  spacing  does  not  appear  to  be 
smooth  in  some  areas.  To  guarantee  that  key  boundary  or  contour  points, 
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Fig.  3.  Contour  and  grid  for  l/|r(z)|. 


Fig.  4.  Virtual  Reality  Modeling  Language  (VRML)  Views  of  l/|r(z)|. 


such  as  those  at  zeros  or  corners,  are  maintained  regardless  of  grid  size, 
we  identify  “fixed  points,”  that  is,  boundary  points  that  must  always  be 
kept.  Grid  lines  are  always  drawn  through  these  points.  In  most  cases, 
the  resulting  discontinuities  in  cell  spacing  do  not  appear  to  affect  the 
quality  of  the  3D  visualizations.  In  Figure  6 we  attract  grid  points  in  an 
equally  spaced  square  grid  to  a circle.  Currently,  we  are  working  on  adap- 
tive movement  based  on  function  curvature,  but  several  problems  must  be 
addressed.  Functional  minimization  routines,  originally  designed  only  for 
n>3  = 0,  are  being  updated  to  accommodate  a likely  nonlinear  dependence 
of  u on  the  spline  coefficients.  Also,  the  specialized  nature  of  many  of  the 
functions  means  that  accessing  and  linking  the  codes  needed  to  compute 
gradient  and  curvature  data  may  not  be  a trivial  task. 

§5.  Conclusions 

We  have  successfully  used  boundary / contour  fitted  grid  generation  to  cre- 
ate over  one  hundred  visualizations  of  complex  mathematical  functions  for 
the  NIST  Digital  Library  of  Mathematical  Functions  Project.  The  grid 
generation  technique  has  helped  us  address  many  of  the  problems  such 
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Fig.  6.  Grid  adapted  to  circular  shape. 


as  poor  clipping,  inaccurate  plots,  and  poor  color  mapping  that  can  come 
with  using  standard  commercial  packages  for  3D  graphics.  We  continue  to 
improve  the  technique  and  have  made  significant  progress  toward  imple- 
menting the  adaptive  movement  of  grid  lines  based  on  external  information 
such  as  function  curvature  and  gradient  data. 

§6.  Disclaimer 

All  references  to  commercial  products  are  provided  only  for  clarification 
of  the  results  presented.  Their  identitication  does  not  imply  recommen- 
dation or  endorsement  by  NIST. 
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